log using replicationlog.log, replace

/*
Replication of "Can't coalesce, can't constrain: Redefining elite influence in non-democracies"
José Kaire
7/16/2022
*/ 

cd "C:\Users\josek\Documents\Research\ECA paper\Replication materials\PSRM replication" 

global format1 graphregion(color(white))  ylabel(, glcolor(gs12) glpattern(dot)) ymtick(##5) xlabel(, grid glcolor(gs12) glpattern(dot) ) xmtick(##5) 


**Figure 2: Simulation results graph 
*Panel A
import delimited "https://github.com/josekaire/ecapaper/raw/main/theta_recovery_random_items.csv", clear
gen id= _n
gen odd = mod(id,2)  
gen marker=.

twoway  (rspike lower90 upper90 id, lcolor(black%90) lwidth(vvthin))  (rspike marker marker  id, lcolor(black%90) lwidth(thin)) (scatter simulatedtheta id, msize(small) msymbol(X) mcolor(black) lcolor(%.40) mlwidth(vthin))  if odd==1&id<400, xsize(6.5) name(g400, replace) legend(order (3 "Simulated dimension, {it:{&theta}{sub:jt}}" 2 "Estimated random items (RI) {it:{&theta}{sub:jt}}") col(1) ring(0) position(11) bmargin(small) region(fcolor(%0) lcolor(white%0))) $format1 xlabel(none) xtitle("") xmtick(none) ytitle("{it:{&theta}{sub:jt}}")

twoway  (rspike lower90 upper90 id, lcolor(black%90) lwidth(vvthin))  (scatter simulatedtheta id, msize(small) msymbol(X) mcolor(black) lcolor(%.40) mlwidth(vthin)) if odd==1&id<800&id>400, xsize(6.5) name(g800, replace) legend(off)  $format1  xlabel(none) xtitle("") xmtick(none) ytitle("{it:{&theta}{sub:jt}}")

twoway  (rspike lower90 upper90 id, lcolor(black%90) lwidth(vvthin))  (scatter simulatedtheta id, msize(small) msymbol(X) mcolor(black) lcolor(%.40) mlwidth(vthin)) if odd==1&id<1200&id>800, xsize(6.5) name(g1200, replace) $format1 legend(off) xlabel(none) xtitle("Observations order by group") xmtick(none) ytitle("{it:{&theta}{sub:jt}}")

graph combine g400 g800 g1200, ycommon col(1) name(panelA, replace) graphregion(color(white)) title("(A): Recovery of latent dimension", col(black) size(medsmall))

*Panel B
import delimited "https://github.com/josekaire/ecapaper/raw/main/rmse.csv", clear 
replace index=4 if parameter=="eta"
encode parameter, gen(param)
gen wrmse=(rmse/weight)*100
twoway (bar wrmse index if parameter=="theta"&model=="m2pl", fcol("black") lcol(black%70) barwidth(.4)) (bar wrmse index if parameter=="theta"&model=="mRI", col("gs10") fcolor(*.7) lcolor(gs10*.7) barwidth(.4)) ///
 (bar wrmse index if parameter=="alpha"&model=="m2pl", fcol("black") lcol(black%70) barwidth(.4)) (bar wrmse index if parameter=="alpha"&model=="mRI", col("gs10") fcolor(*.7) lcolor(gs10*.7) barwidth(.4)) ///
  (bar wrmse index if parameter=="beta"&model=="m2pl", fcol("black") lcol(black%70) barwidth(.4)) (bar wrmse index if parameter=="beta"&model=="mRI", col("gs10") fcolor(*.7) lcolor(gs10*.7) barwidth(.4)) ///
	(bar wrmse index if parameter=="eta"&model=="m2pl", fcol("black") lcol(black%70) barwidth(.4)) (bar wrmse index if parameter=="eta"&model=="mRI", col("gs10") fcolor(*.7) lcolor(gs10*.7) barwidth(.4)), legend(order (1 "2PL model" 2 "RI model") col(1) ring(0) position(11) bmargin(small) region(fcolor(%0) lcolor(white%0))) $format1 title("(B): Error in 2PL and RI models", col(black) size(medsmall))xmtick(none) ymtick(none)   yscale(range(90 150)) ytitle("Root mean square error (%)")   xlabel( 1 "{it:{&theta}{sub:jt}}" 2 "{it:{&alpha}{sub:k}}" 3 "{it:{&beta}{sub:k}}" 4 "{it:{&eta}{sub:jtk}}") xtitle("") xscale(range(.5 4.5)) name(rmse, replace)

*Panel C 
import delimited "https://github.com/josekaire/ecapaper/raw/main/eta_recovery.csv", clear
gen marker=.
sort simulated 
gen invsimulated=invlogit(simulated)
gen invm1=invlogit(pleta)
gen invm2=invlogit(ranitemeta)

twoway (scatter  invm1 simulatedeta, msymbol(o) mlcolor(black%0) mcolor(black%70))(scatter  invm2 simulatedeta, msize(small) mcolor(gs16%20) mlcolor(black%30) mlwidth(vvthin) msymbol(d)) (scatter marker simulatedeta, mcolor(gs13%50) mlcolor(black%70) mlwidth(vvthin) msymbol(d))(line invsimulated simulated, lcolor(black%70)) if simulated>-10&simulated<10, title("(C): Overall model performance", col(black) size(medsmall))ytitle("Probability of a positive response") xtitle("True {it:{&eta}{sub:jtk}}") $format1 legend(size(small) order (4 "True probability" 1 "2PL model" 3 "RI model" ) col(1) ring(0) position(11) bmargin(small) region(fcolor(%0) lcolor(white%0))) xscale(range(-8 8)) xmtick(none) name(etarecover, replace)

graph combine rmse etarecover, name(BandC, replace) col(1) graphregion(color(white))
graph combine panelA BandC, graphregion(color(white))
graph export figure2.pdf, replace

*Figure 3: Correlation of power-sharing and cohesion 
use  "https://github.com/josekaire/ecapaper/raw/main/Replication%20data.dta", clear 
twoway (scatter dps bpscondensed , mcolor("gs3") mfcolor(%5) mlcolor(%0)) (qfit dps bps, lcolor("black")) if bps>0, graphregion(color(white))  ylabel(, glcolor(gs12) glpattern(dot) nogextend) ymtick(##5) xlabel(, grid glcolor(gs12) glpattern(dot) nogextend) xmtick(##2)  ytitle("Elite power") xtitle("Elite cohesion") legend(off) title("(A): Correlation of elements of collective action", color("black") size("medium")) name(pooled, replace) 

preserve
import delimited "https://github.com/josekaire/ecapaper/raw/main/power%20and%20cohesion%20correlation%20by%20country.csv", varnames(1) clear 
destring rhos uniques sds, replace ignore(`"NA"')
drop if rhos==.
gen weighted=rhos*(uniques)
sort weighted 
gen id=_n
gen cumfreq = id/_N
twoway spike weighted cumfreq, lcolor("black")  graphregion(color(white))  ylabel(, glcolor(gs12) glpattern(dot) nogextend) ymtick(##5) xlabel(, grid glcolor(gs12) glpattern(dot) nogextend) xmtick(##2)  ytitle("Correlation power and cohesion") xtitle("Cumulative frequency by regime") title("(B): Correlation by regime", color("black") size("medium")) name(byregime, replace)
restore 

graph combine pooled byregime ,  graphregion(color(white)) name(eca, replace) col(2) xsize(7) 	
graph export figure3.pdf, replace

*Figure 4: ECA within reigme types 
**Variations within type 
reg eca  if country=="Mexico"&year>1957
margins,  post
estimates store mex
reg eca  if country=="Malaysia"&year>1957
margins,  post
estimates store mal
reg eca  if country=="Indonesia"&year>1957
margins,  post
estimates store indo
reg eca  if country=="Morocco"&year>1955
margins,  post
estimates store mor
reg eca  if country=="Oman"&year>1955
margins,  post
estimates store oman
reg eca  if country=="Nepal"&year>1955
margins,  post
estimates store nep

coefplot (mex,label("Mexico") col("black"))  (mal,label("Malaysia") col("gs5") ) (indo,label("Indonesia") col("gs10")), vertical recast(bar) barwidth(0.1)  fcolor(*.5)  ciopts(recast(rcap)col("black")) citop title("Party regimes", col("black")) graphregion(color(white))   xlabel( .75 "Mexico" 1 "Malaysia" 1.25 "Indonesia" ,grid glcolor(gs12) glpattern(dot))  grid(between glcolor(orange) glpattern(dash)) legend(off) name(parties, replace)
  
coefplot (mor,label("Morocco") col("black"))  (nep,label("Nepal") col("gs5") ) (oman,label("Oman") col("gs10")), vertical recast(bar) barwidth(0.1)  fcolor(*.5)  ciopts(recast(rcap)col("black")) citop title("Monarchies", col("black")) graphregion(color(white))   xlabel( .75 "Morocco" 1 "Nepal" 1.25 "Oman" ,grid glcolor(gs12) glpattern(dot))  grid(between glcolor(orange) glpattern(dash)) legend(off)   name(monarchies, replace)  

graph combine parties monarchies, ycommon xsize(7.5)  graphregion(color(white))
graph export figure4.pdf, replace

*Figure 5: Mexico Zimbabwe comparison  
qui twoway (tsline eca if ccode==70, lcolor(black)) ///
 (tsline gl if ccode==70, lcolor(gray)) ///
 (tsline gwf_ps if ccode==70, lcolor(gray) lpattern(dash)), ///
 graphregion(color(white)) ylabel(, glcolor(gs10) glpattern(dot)) xlabel(#8, grid glcolor(gs10) glpattern(dot)) xmtick(##5) xtitle("Year") ytitle("Elite influence") title("Mexico", color("black"))  name(mex, replace)  legend(ring(0) position(5)  label(1 "ECA") label(2 "GS") label(3 "GWF") col(1) region(style(none)))    
 
qui twoway (tsline eca if ccode==552, lcolor(black)) ///
 (tsline gl if ccode==552, lcolor(gray) ) ///
 (tsline gwf_ps if ccode==552, lcolor(gray) lpattern(dash)) if year>1979, graphregion(color(white))  ylabel(, glcolor(gs10) glpattern(dot)) xlabel(#8, grid glcolor(gs10) glpattern(dot)) xmtick(##5) xtitle("Year") ytitle("Elite influence")  title("Zimbabwe", color("black")) name(zimb, replace) legend(ring(0) position(5) label(1 "ECA") label(2 "GS") label(3 "GWF") col(1) region(style(none)))

graph combine mex zimb, xsize(7) ycommon  graphregion(color(white))
graph export figure5.pdf, replace

*Table 4: Replication of Gehlbach and Keefer 
qui reg fdi4_gdp gov1_yrs tenure stabs fuelex_gdp oresex_gdp frac_ethn frac_relig frac_ling pop_yng_pct pop_tot pop_ru_pct land_km gdppc_ppp_2005_us if eca!=., cluster(ifs)  


qui reg fdi4_gdp eca tenure stabs fuelex_gdp oresex_gdp frac_ethn frac_relig frac_ling pop_yng_pct pop_tot pop_ru_pct land_km gdppc_ppp_2005_us  if insample!=., cluster(ifs) 
estimates store m3 

qui reg fdi4_gdp gov1_yrs_std tenure stabs fuelex_gdp oresex_gdp frac_ethn frac_relig frac_ling pop_yng_pct pop_tot pop_ru_pct land_km gdppc_ppp_2005_us if eca!=., cluster(ifs) 
estimates store m2

qui reg fdi4_gdp gov1_yrs_std tenure stabs fuelex_gdp oresex_gdp frac_ethn frac_relig frac_ling pop_yng_pct pop_tot pop_ru_pct land_km gdppc_ppp_2005_us, cluster(ifs) 
estimates store m1

estout  m1 m2 m3,  label  cells(b(star fmt(2))& se(par fmt(2))) s(N r2)

*Figure 6: Human rights and elite influence 
*Rescale Fariss measure 
sum hrs
gen HRS=hrs+abs(r(min))
sum HRS
replace HRS=HRS/r(max)
sum HRS
drop if sid2==.
xtset sid2 year 
foreach var of varlist eca gl_ps gwf_ps {
	forv j=4(1)10{
		forv i=1(1)500{
				*Run regression 
				gen rng`i'=uniform()
				qui xtreg HRS `var' gwf_pa gwf_mil gwf_mon if rng`i'<(`j'/10)
				*Save coefficient/std error
			   
			   if `i'==1 matrix `var'_`j'=_b[`var']/_se[`var'] 
			   else matrix `var'_`j'=(`var'_`j' \ _b[`var']/_se[`var'] )
			   
			   drop rng`i'
			   
		}
	display `j'
	svmat `var'_`j'
	}
}
preserve
*reshape data
keep eca_41-gwf_ps_101
gen id=_n
reshape long eca_ gl_ps_ gwf_ps_ , i(id) j(level)
drop if eca==.
gen jitter=rnormal(0, .0001)
replace gwf_ps=gwf_ps+jitter if level==101
replace gl_ps=gl_ps+jitter if level==101
replace eca=eca+jitter if level==101
replace level=level-1
*Produce box plots
graph box eca, over(level)  box( 1 ,col(black%50) lcolor(black%70)) marker(1, mfcolor(black%50) mlcolor(black%0)) ylabel(, glcolor(gs12) glpattern(dot)) ymtick(##5)  graphregion(color(white)) title("ECA measure", color("black")) ytitle("t-value") yline(1.9, lc(black%80)) text( 2.2 100  "{it:p}<.05") name(eca, replace) b1title("Percent of sample used")
graph box gl_ps, over(level)  box( 1 ,col(black%50) lcolor(black%70)) marker(1, mfcolor(black%50) mlcolor(black%0)) ylabel(, glcolor(gs12) glpattern(dot)) ymtick(##5)  graphregion(color(white)) title("GS measure", color("black")) ytitle("t-value") yline(1.9, lc(black%80))  name(gs, replace) b1title("Percent of sample used")
graph box gwf_ps, over(level)  box( 1 ,col(black%50) lcolor(black%70)) marker(1, mfcolor(black%50) mlcolor(black%0)) ylabel(, glcolor(gs12) glpattern(dot)) ymtick(##5)  graphregion(color(white)) title("GWF measure", color("black")) ytitle("t-value") yline(1.9, lc(black%80)) name(gwf, replace) b1title("Percent of sample used")
graph combine eca gs gwf, ycommon col(3) xsize(7) graphregion(color(white)) name(graph1, replace) 
graph save graph1, replace
*
restore 

*AIC comparison	*
xtset sid2 year 

foreach var of varlist eca concentrationofpower latent_personalism {
		*Run regression 
		qui xtreg hrs `var' gwf_pa gwf_mil gwf_mon if eca!=.&gl_ps!=.&gwf_ps!=., mle
		qui estat ic
		*Save coefficient/std error
		mat s=r(S)	   
		if `var'==eca matrix AIC=s[1,5]
	    else matrix AIC=(AIC \ s[1,5] )	   	
	}

svmat AIC
sum AIC1, meanonly
replace AIC1=AIC1-r(min)+1
gen id=_n if _n<4
graph bar (mean) AIC1, over(id, relabel(1 "ECA" 2 "GS" 3 "GWF")) box( 1 ,col(black%50) lcolor(black%70)) marker(1, mfcolor(black%50) mlcolor(black%0)) ylabel(, glcolor(gs12) glpattern(dot)) ymtick(##5)  graphregion(color(white)) ytitle("Relative Error (AIC)") title("(A) Predictive power of elite influence measures", col("black") size(medium)) name(aic, replace)

graph combine aic  graph1 , col(1) graphregion(color(white))  
graph export figure6.pdf, replace

**Online appendix 
*Item parameters, Figure 1A
import delimited "https://github.com/josekaire/ecapaper/raw/main/alpham1plot.csv", clear
rename v1 id
twoway  (rspike lowerhpdi50 higherhpdi50 id, lcolor(black%60)) (rspike lowerhpdi70 higherhpdi70 id, lcolor(black%30) lwidth(medthin)) (rspike lowerhpdi90 higherhpdi90 id, lcolor(black%15)  lwidth(thin))  (scatter pointestimate id, msymbol(o) msize(medsmall) mcolor("gs11")) (scatter alpha id, msize(small) msymbol(X) mcolor(black)), title("(A): 2PL {it:{&alpha}{sub:j}} recovery", col(black) size(medsmall)) $format1  xtitle("Item ID") ytitle("Discrimination parameter, {it:{&alpha}{sub:j}}") legend(size(small) order (4 "Estimated {it:{&alpha}{sub:j}}" 5 "True {it:{&alpha}{sub:j}}") col(1) ring(0) position(11) bmargin(small) region(fcolor(%0) lcolor(white%0)))  xlabel(#10) xmtick(##1) name(alpham1, replace)

import delimited "https://github.com/josekaire/ecapaper/raw/main/alpham2plot.csv", clear
rename v1 id
twoway  (rspike lowerhpdi50 higherhpdi50 id, lcolor(black%60)) (rspike lowerhpdi70 higherhpdi70 id, lcolor(black%30) lwidth(medthin)) (rspike lowerhpdi90 higherhpdi90 id, lcolor(black%15)  lwidth(thin))  (scatter pointestimate id, msymbol(o) msize(medsmall) mcolor("gs11")) (scatter alpha id, msize(small) msymbol(X) mcolor(black)), title("(B): RI {it:{&alpha}{sub:j}} recovery", col(black) size(medsmall)) $format1  xtitle("Item ID") ytitle("Discrimination parameter, {it:{&alpha}{sub:j}}") legend(size(small) order (4 "Estimated {it:{&alpha}{sub:j}}" 5 "True {it:{&alpha}{sub:j}}") col(1) ring(0) position(11) bmargin(small) region(fcolor(%0) lcolor(white%0))) xlabel(#10) xmtick(##1) name(alpham2, replace)

import delimited "https://github.com/josekaire/ecapaper/raw/main/beta_m1plot.csv", clear
rename v1 id
twoway  (rspike lowerhpdi50 higherhpdi50 id, lcolor(black%60)) (rspike lowerhpdi70 higherhpdi70 id, lcolor(black%30) lwidth(medthin)) (rspike lowerhpdi90 higherhpdi90 id, lcolor(black%15)  lwidth(thin))  (scatter pointestimate id, msymbol(o) msize(medsmall) mcolor("gs11")) (scatter beta id, msize(small) msymbol(X) mcolor(black)), title("(C): 2PL {it:{&beta}{sub:j}} recovery", col(black) size(medsmall)) $format1  xtitle("Item ID") ytitle("Difficulty parameter, {it:{&beta}{sub:j}}") legend(size(small) order (4 "Estimated {it:{&beta}{sub:j}}" 5 "True {it:{&beta}{sub:j}}") col(1) ring(0) position(11) bmargin(small) region(fcolor(%0) lcolor(white%0))) xlabel(#10) xmtick(##1) name(betam1, replace)

import delimited "https://github.com/josekaire/ecapaper/raw/main/beta_m2plot.csv", clear
rename v1 id
twoway  (rspike lowerhpdi50 higherhpdi50 id, lcolor(black%60)) (rspike lowerhpdi70 higherhpdi70 id, lcolor(black%30) lwidth(medthin)) (rspike lowerhpdi90 higherhpdi90 id, lcolor(black%15)  lwidth(thin))  (scatter pointestimate id, msymbol(o) msize(medsmall) mcolor("gs11")) (scatter beta id, msize(small) msymbol(X) mcolor(black)), title("(D): RI {it:{&beta}{sub:j}} recovery", col(black) size(med)) $format1  xtitle("Item ID") ytitle("Discrimination parameter, {it:{&beta}{sub:j}}") legend(size(small) order (4 "Estimated {it:{&beta}{sub:j}}" 5 "True {it:{&beta}{sub:j}}") col(1) ring(0) position(11) bmargin(small) region(fcolor(%0) lcolor(white%0)))  xlabel(#10) xmtick(##1) name(betam2, replace)

graph combine alpham1 alpham2, ycommon name(alpha, replace) graphregion(color(white))
graph combine betam1 betam2, ycommon name(beta, replace) graphregion(color(white))

graph combine alpha beta, graphregion(color(white)) col(1)
graph export figure7.pdf, replace
*Item random effects, Figure 2A
import delimited "https://github.com/josekaire/ecapaper/raw/main/realdata_looic.csv", encoding(UTF-8)  clear
encode model, gen(mri)
encode parameter, gen(power)

qui reg looic i.mri##i.power 
 qui margins, at(mri=(1 ) power=(1 2)) post
qui estimates store m2PL 
qui reg looic i.mri##i.power 
qui margins, at(mri=(2) power=(1 2)) post
estimates store mRI 
 
coefplot (m2PL, fcolor(black%80) ) (mRI, fcolor(gray%80)) ,  coeflabels(1._at ="Power-sharing" 2._at ="Elite cohesion")  vertical recast(bar) barwidth(.3)   fcolor(.5) lcolor(%0) legend(size(small) order (2 "2PL model" 4 "RI model") col(1) ring(0) position(11) bmargin(small) region(fcolor(%0) lcolor(white%0))) yscale(range(3600 4800)) ylabel(3800(200)4600)  graphregion(color(white))  ylabel(, glcolor(gs12) glpattern(dot)) ymtick(##0)  xlabel(, grid glcolor(gs12) glpattern(dot) ) xlabel(, grid glcolor(gs12) glpattern(dot) ) 
graph export figure8.pdf, replace

log close
graphlog using replicationlog.log, lspacing(1) replace
